T/V- 76> 


NASA Technical Memorandum 86412 nasa-tm-86412 19850017209 


CRYSTALLINITY DETERMINATION BY CURVEFIT PROCEDURE FOR A SEMI- 
CRYSTALLINE POLYMER 


N. T. WAKELYN 


MAY 1985 



NASA 

National Aeronautics and 
Space Administration 

Langley Research Center 

Hampton, Virginia 23665 




SUMMARY 


Wide angle x-ray scattering (WAXS) data from poly(etheretherketone) 
(PEEK) has been resolved into a crystalline contribution represented as 4 
reflections and an amorphous contribution represented as a broad, smoothly 
varying curve, both contributions occurring in the 2e range: 15-31 

degrees. In this resolution the crystalline scatter is described as a 
linear combination of Cauchy and Gaussian functions while that of the 
amorphous halo is expressed as a cubic polynomial. Statistical analysis of 
the measured scattered intensity from an amorphous specimen with that 
calculated from the cubic polynomial, as a function of the combination 
parameter (fraction of Cauchy and Gaussian functions), suggests that the 
crystalline fraction of the polymer specimen studied is about 0.39. A 
listing of the FORTRAN IV program used in the resolution is provided in the 
Appendix. 


INTRODUCTION 


Composites fabricated from thermoplastic resins have some inherent 
differences in physical properties from those of the current thermoset resin 
composites and some formulations show considerable promise for structural 
applications. A partial listing of some of the potential advantages of such 
materials should include ease of formability, damage repair, and bonding as 
well as good environmental resistance and re-processing capability. 1-3 
The concept of a thermoplastic implies the possibility of some degree of 
crystallization and with this possibility there is the promise of material 
property enhancement. Indeed, increased toughness, reduced moisture 






degradation, and the chemical resistance required for use in aerospace 
structures have been associated with crystallization in thermoplastic resin 
matrices. 4 

Poly(etheretherketone) (PEEK) is one of the thermoplastic matrix resins 
which can be fabricated as a semi -crystal line matrix in combination with a 
carbonaceous fiber reinforcement. This material has been characterized as 
having mechanical properties superior to those of currently used epoxies 
when tested under wet conditions at elevated temperatures. 4 Physical 
property data and details of the synthesis of PEEK and other polyaryl ether- 
ketones have been reported. 5 * 6 

In order to understand the relationship between a desired physical 
property and an induced crystalline order in a matrix material a method for 
estimation of the degree of this order is necessary. Added complexity in a 
polymeric structure, e.g. changing from a one to a two phase material, is 
sure to have a far from straightforward effect on desired physical 
properties. It is well within reason to expect a small amount of added 
crystalline order to produce an increase in toughness, while a greater 
increase yields a brittle material. Hopefully a method for measuring 
crystalline order in a polymer matrix would also be useful in the presence 
of a carbonaceous reinforcement. 

Wide angle x-ray scattering (WAXS) has been used to investigate the 
nature of the crystalline structure appropriate for the poly( aryl ether) and 
poly(aryletherketone) polymers. 7-9 The basis for the interpretation of 
the WAXS data for PEEK and the other poly(aryletherketone) polymers has been 
the study of poly( p-phenyl ene oxide) (PPO) which is reported 7 to crystallize 
in the space group Pbcn (space group No. 60, International Tables for X-ray 
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Crystallography). While PEEK is a more complicated polymer than PPO and 
requires at least six aryl units aligned along the c-axis of its unit cell 
for its description, 6 * 8 the symmetry elements of Pbcn are appropriate for 
PEEK, as well as for PPO, and the major peaks in the WAXS data for this 
material can be indexed based upon the assumption of a structural analogy 
between the two materials. 10 The four major crystalline reflections of 
PEEK, occurring at -18.7, 20.7, 22.6, and 28.7° (2a), have, therefore, 
been indexed as 110, 113, 200, and 213, respectively, in strict analogy with 
the indexing of reference 9. 

The purpose of this report is to describe an analytical procedure for 
the estimation of degree of crystalline order in a thermoplastic polymer 
matrix. Semicrystal 1 1 ne PEEK, a material of possible future use in 
aerospace applications, has been used to develop and illustrate this 
procedure. Intensity versus diffraction angle data taken from a flat powder 
specimen in a WAXS experiment has been used for the analysis. In order to 
build in a certain realism and thus to facilitate extention to real 
specimens of composites reflection geometry has been used to obtain the 
data. While the computer program developed for the analysis and presented 
in the Appendix is specialized for PEEK it is readily convertible to other 
formulations for which the major x-ray diffraction reflections are known. 


ANALYTICAL 

Various methods for the determination of degree of crystalline order in 
polymeric materials using WAXS data have been employed, some of which have 
achieved an almost standard practice status. A sampling of these methods. 
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in order of increasing ease of application, vary from a difficult to use but 
sophisticated absolute method, 11 through a more practical though less 
rigorous absolute procedure, 12 through a well -regarded index 
procedure, 13 ’ 14 to a surprisingly useful 2-point index method. 15 All of 
these would be confounded by the presence of the carbonaceous reinforcement 
and involve factors such as a requirement for multiple samples, transmission 
geometry, interpolation between standards, and for some, even an element of 
subjectivity. For those polymers with known crystal structure, however, a 
crystallinity determination by curvefit procedure 16 is feasible and in 
theory, at least, is unaffected by the previously mentioned factors and may 
even be practicable in the presence of the carbonaceous reinforcement. 

The basis of the curvefit procedure 16 is the expression of the measured 
x-ray intensity data versus diffractometer angle (26) for the polymeric 
specimen as a sum of the x-ray scatter from the known crystalline peaks plus 
that characteristic of the broad, smoothly varying paracrystalline or 
amorphous scatter. The crystalline peaks are described as a linear 
combination of Cauchy and Gaussian expressions with the linear combination 
parameter varying from 0 for a pure Cauchy function to 1 for a pure Gaussian 
function with intermediate values representing a combination of the 2 
functions. The amorphous background scatter is expressed as a cubic 
polynomi al . 

Each crystalline peak is defined by parameters representing the peak 
height, peak position, and width of the peak at half maximum, while the 
cubic polynomial is defined by 4 parameters. Thus, the expression for the 
summation of the crystalline peaks plus that of the background curve for the 
4 reflections representative of PEEK in the 20 range of interest contains 16 
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independent variables to be determined by the curvefit procedure. These 
independent variables were determined by solving the system of nonlinear 
equations, i.e. the summation of the expressions for the crystalline 
reflections and background curve as a function of the diffraction angle 
minus the measured intensity value at that angle, by using Levenberg- 
Marquardt and Gauss algorithms for the nonlinear least squares approxima- 
tions . 17 The computer program, which solves for these variables and which 
contains the Langley computing center's library routine, MARQ , 18 is present- 
ed in the Appendix for a real case. Convergence of the modified Marquardt 
algorithm is satisfied when the difference of the residual sum of squares 
estimates on two successive iterations divided by the residual at the first 
of these iterations is equal to as less than an input parameter based upon 
the relative accuracy of the equations . 18 

Once the parameters have been obtained, an integration or comparable 
mathematical procedure may be used to extract the relative areas under the 
various curves. An obvious first approximation to the polymeric crystalline 
fraction would then be the summation of the areas under the crystalline 
reflections divided by the area under all curves. There is apparently some 
controversy about this, however, since some workers 12 report different 
constants of proportionality for relating the crystalline and amorphous 
fractions to the respective integral intensities (areas). Other workers , 19 
base their data reduction upon the "Law of Conservation of Intensity ", 20 in 
effect, choosing the more simple procedure. Since the present goal is the 
description of an analytical procedure as applied to a specific 
thermoplastic material and since only one crystalline specimen was required 
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to develop and illustrate this procedure, the more simple expression for 
crystalline fraction is presented. 


EXPERIMENTAL 


WAXS data has been obtained on annealed and melt-quenched PEEK 
material. CuKa radiation was used with an automated powder x-ray diffracto- 
meter equipped with a curved crystal, graphite monochromator. With the gen- 
erator operated at 45 kV and 40 mA, the intensity of 1 s counts taken every 
0.01 degrees (20) was recorded on hard disk. The data in the angular range 
15.00-30.83°, an adequate range for evaluation of crystalline fraction, 12 
was used for the present analysis. 

Commercially synthesized polymeric material, received in the form of 
compression mouldings, was melted in a hot press and immediately quenched in 
ice water to produce a light brown, transparent film which had a density of 
1.267 g cm -3 as determined by immersion in a toluene-carbon tetrachloride 
density gradient column. This film was cooled to dry ice temperature, 
ground in a rotary mill, and sieved using a 60 mesh screen. Part of the 
material which passed through the screen was used as is (the melt-quenched 
specimen) and part was annealed in a vacuum oven at 267°C for 72 h to induce 
crystallinity. The density of this annealed material (the crystalline 
specimen) was determined to be 1.295 g cm -3 . 

RESULTS AND DISCUSSION 


The mathematical expression used to describe the crystalline 
reflections contains a linear combination parameter which controls their 
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geometrical shape. When this combination parameter has the value 0, the 
reflections are characterized as pure Cauchy functions, and when its value 
is 1 they are represented as pure Gaussians. Intermediate values of the 
parameter produce composite representations of the reflections. Figure 1 
presents a series of 10 plots (Figures 1(a) - l(k)) depicting the changes 
incurred in the crystalline peaks and the background curve as the linear 
combination parameter is varied from 1.0 to 0.0 in 0.1 increments. They are 
compared on each plot with the actual measured data of semi -crystal 1 i ne PEEK 
appearing at the top of the figure as a curve composed of a series of points 
over the angular range 15.00-30.83° (28). The indices describing the 
crystalline reflections are included in Figure 1(a). 

As the description of the crystalline peaks changes from pure Gaussian 
(Fig. 1(a)) to pure Cauchy (Fig. l(k)) tailing at the base of the peaks 
increases uniformly. This increased tailing produces a directly related 
increase in the area of the crystalline peaks, which increase at the expense 
of the background area. The fractional crystalline area, which is the sum 
of areas under the 4 crystalline peaks divided by this sum plus the area 
under the background curve, varies with the increased tailing from 0.2502 to 
0.3916 in arbitrary intensity-degree units. This fractional crystalline 
area is related to a degree of crystalline order with the assumption of 
equivalent scatter from equivalent phases. In all cases there is 
insignificant difference between the calculated area under the measured 
curve and the calculated summation of all curvefit areas. 

Since the summation of crystalline areas plus the area under the 
background curve is essentially constant, and since the summation of the 
crystalline areas changes with the linear combination parameter, the shape 
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of the background curve (as well as its area) also changes with this 
parameter. A comparison of this background curve with that measured for 
amorphous PEEK, normalized for equal area, is presented in Figure 2, as a 
function of the linear combination parameter. There is a progressive 
observable difference in the relationship between the 2 compared curves 
presented in Figures 2(a) and 2(b); however, a statistical analysis was 
performed to minimize the necessity for subjective judgment as to the 
quality of a representation. The simple correlation coefficients between 
the set of points describing the background curve as a function of the 
linear combination parameter and the measured data for an amorphous 
specimen, each having 1584 observations, were determined and are presented 
in Table I. 

The variation of the amorphous area, fractional crystalline area, and 
correlation coefficient as the form of the crystalline peak shape is 
systematically varied from that of a pure Cauchy function to that of a pure 
Gaussian function is presented in Table I. There is displayed a systematic 
increasing trend in the correlation coefficient with increasing Cauchy 
character suggesting that for the present WAXS experiment with 
semi -crystal 1 1 ne PEEK powder the pure Cauchy description of the crystalline 
peak shape is preferred. Others, however, have reported a linear 
combination parameter of 0.5 to yield the best fit for many synthetic and 
natural fibers, 16 while another worker analyzed nylon data using pure Cauchy 
expressions. 21 Concomitant with the present trend is the determination that 
the crystalline fraction of the annealed PEEK is 0.39. 
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CONCLUDING REMARKS 


A curvefit procedure, based upon the modified Marquardt algorithm, has 
been developed and illustrated with semi -crystal line PEEK. Intensity versus 
diffraction angle data taken from flat powder specimens (crystalline and 
amorphous) in reflection geometry has been used for the analysis. While the 
method was designed to be based upon one semi -crystal line sample only, an 
amorphous specimen was required to discriminate between the various 
solutions. Extention of the analytical procedure (and computer program) to 
other systems, such as PEEK composites containing carbonaceous fiber, should 
be straightforward. 
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Table I. Results of Curvefitting Procedure and Statistical Analysis as a 

Function of the Linear Combination Parameter for Semi -crystal line PEEK 
Polymer Matrix. 




FRACTIONAL 


COMBINATION 

AMORPHOUS 

CRYSTALLINE 

CORRELATION 

PARAMETER* 

AREA** 

AREA 

COEFFICIENT*** 

0.0 

23441 

.3916 

.9705 

0.1 

24079 

.3751 

.9701 

0.2 

24690 

.3592 

.9698 

0.3 

25276 

.3440 

.9694 

0.4 

25842 

.3293 

.9692 

0.5 

26388 

.3151 

.9690 

0.6 

26917 

.3014 

.9687 

0.7 

27431 

.2881 

.9685 

0.8 

27930 

.2751 

.9683 

0.9 

28416 

.2625 

.9681 

1.0 

28891 

.2502 

.9679 


*A linear combination parameter of 0.0 implies a pure Cauchy function 
while that of 1.0 implies a pure Gaussian expression. 

**Area under amorphous halo in intensity - degree units relative to an 
area under the measured intensity envelope of 38530. 

***A measure of how well the shape of the derived cubic polynomial, 

describing the amorphous halo, is correlated to the measured shape of 
an amorphous PEEK specimen. 
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Figure 1(a). - X-ray diffractogram of semi-crystalline PEEK compared to resolved 
crystalline scatter represented as a linear combination of Cauchy and Gaussian 
functions and the amorphous halo represented by a cubic polynomial. 
Combination parameter = 1.0. Fractional crystalline area = 0.2502. 



Figure 1(b). - Combination parameter = 0.9. Fractional crystalline area = 0.2625. 
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Figure 1(c).- Combination parameter = 0.8. Fractional crystalline area = 0.2751. 



Figure 1(d). - Combination parameter = 0.7. Fractional crystalline area = 0.2881. 
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Figure 1(g). - Combination parameter = 0.4. Fractional crystalline area = 0.3283. 



Figure 1(h). - Combination parameter = 0.3. Fractional crystalline area = 0.3440. 
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Figure l(i). - Combination parameter = 0.2. Fractional crystalline area = 0.3592. 



Figure l(j). - Combination parameter = 0. 1. Fractional crystalline area = 0.3751. 
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Figure 2(a). - Amorphous halo as represented by a cubic polynomial compared 
with data from amorphous PEEK, normalized for equal area, for combination 
parameters varying from 0.5 to 1.0. 
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Figure 2(b). - Amorphous halo as represented by a cubic polynomial compared 
with data from amorphous PEEK, normalized for equal area, for combination 
parameters varying from 0. 0 to 0. 4. 
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APPENDIX I. A FORTRAN IV COMPUTER PROGRAM FOR THE 
RESOLUTION OF POLYMER MATRIX CRYSTALLINE SCATTER 


The code listed in this Appendix represents that used to calculate one of. 
the cases presented in this report. A dummy call to a PLOT routine is used 
since these routines are, in general, arbitrary and thus not transportable. 
Subroutine MARQ, which solves a system of N nonlinear equations in N 
unknowns for the real roots by the modified Marquardt algorithm, is a 
Langley Research Center library routine. (Mathematical and Statistical 
Software at Langley, Central Scientific Computing Complex Document N-3, 
April 1984). 
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PROGRAM MARFIT( INPUT, OUTPUT, TAPE5=INPUT,TAPE6=0UTPUT) 
EXTERNAL F 
COMMON TT, XRI , G 

DIMENSION X( 16) , FX(1584), WK(33800), TT(1584), XRI(1584) 
DIMENSION IRAY( 6) 


DIMENSION IN( 24) 


DATA I RAY/ 6 * (-0)/ 
DATA M/1584/, N/16/, 

C/1.0/, EPS/1 .0E-08/ , 

NSIG/5/, ITMAX/100/ 

DATA 

X/ 2766., 

18.71, 

1.0, 

1 

952. , 

20.67, 

1.0, 

2 

1570., 

22.56, 

1.0, 

3 

1008. , 

28.65, 

1.0, 

5 

-9500. , 

1300. , 

-50. , 

6 

COMMENT. 

0.600/ 

INITIALIZE G AT 

THIS TIME 



G = 0.5 

COMMENT: OBTAIN X-RAY INTENSITIES FROM ARBITRARY FILE 
DO 5 1=1, M 
READ( 5 , 3) IN 

3 F0RMAT(8(A1,A5,A4)) 

C IF( IN(2) . EQ. 5H10.00) GO TO 6 
IF(IN(2) .EQ. 5H14.98) GO TO 6 

5 CONTINUE 

6 PRINT4 , IN(2) 

4 FORMAT( IX, A5) 

11 = 1 

12 = II + 7 

DO 25 13 = 1,6 

DO 50 14 = 1,31 

READ(5, 7) (XRI(I),I=I1,I2) 

7 FORMAT( I0X,8(F7.0)) 

11 = II +8 

12 = 12 + 8 
50 CONTINUE 

12 = 12 - 6 

READ( 5,8) (XRI(I) ,1=11,12) 

8 FORMAT( I0X,2(F7.0) ) 

11 = II + 2 

12 = II + 7 

25 CONTINUE 

DO 26 J=1501, 1584,8 
KEND=J+7 

IF( J .EQ. 1581)KEND=1584 
READ( 5,7) (XRI(I) ,I=J,KEND) 

26 CONTINUE 

PRINT*, KEND 

COMMENT: WE NOW HAVE 1584 XRI'S 
IRAY(4) = 0 

CALL SYSTEMC (115.IRAY) 

PRINT1 

1 FORMAT( 1H1 , I8X, 17HEXPERIMENTAL DATA) 

PRINT2 

2 F0RMAT( 1H0 , 13X, 9H 2-THETA , 10X, 9HINTENSITY, //) 
COMMENT: GENERATE 2-THETA'S CORRESPONDING TO XRI'S 



T = 14.99 
DO 100 1=1, M 
T = T + 0.01 
TT( I) = T 
100 CONTINUE 
105 CONTINUE 

COMMENT: MARQ USES GAUSS AND LEVENBERG-MARQUARDT ALGORITHMS 
C TO RESOLVE WAXS INTO CRYSTALLINE AND PARACRYSTALLINE CONTRIBUTIONS 
CALL MARQ(M,N,X,C,EPS,NSIG,ITMAX,FX,F,WK,IERR) 

PRINT150.G 

150 FORMAT( 1H0 , 3HG =,G22.14) 

PRINT900 , IERR 

IF( IERR .NE. 0) GO TO 10 

C WRITE(6,901) (X( I) , 1=1 ,N) , (FX( I) , 1=1 ,M) , ITMAX 
PRINT200 

200 FORMAT( 1H1 , 25X, 15HS0LUTI0N VECTOR,//) 

PRINT250 

250 FORMAT( 1H0 , 5X, 1 1HPEAK HEIGHT, 10X, 13HPEAK POSITION, 7X, 17HWIDTH AT 1 
1/2 MAX.,//) 

NX = N - 4 

PRINT300, (X(I) ,1=1, NX) 

300 FORMAT ( 1H ,3G22.14) 

PRINT325 

325 FORMAT( 1H0 , 35X, 23HPOLYNOMIAL COEFFICIENTS) 

NP = NX+1 

PRINT350, (X(I), I=NP,N) 

350 FORMAT( 1H0,4G22. 14) 

PRINT400 

400 FORMAT( 1H1 , 30X, 19HVALUES OF EQUATIONS,////) 

PRINT300 , ( FX( I) , 1=1 ,M) 

PRINT500, ITMAX 

500 FORMAT( 23H0NUMBER OF ITERATIONS =,I4) 

CALL AREA(X,M,N) 

CALL PLOT(X,M,N) 

GO TO 20 

10 WRITE( 6 ,902) (X(I),I=1,N) 

PRINT500, ITMAX 
20 STOP 

900 FORMAT( 13H0ERR0R CODE =,I4) 

901 FORMAT( 18HOSOLUTION VECTOR = , 3G22. 14/22H0VALUES OF EQUATIONS =, 

1 4G22 . 14/23H0NUMBER OF ITERATIONS =,I3) 

902 FORMAT( 33H0LAST APPROXIMATION TO SOLUTION =,3G22.14) 

END 

SUBROUTINE PLOT(X,M,N) 

RETURN 

END 1 

SUBROUTINE F(X,M,N,FX) 

COMMON TT, XRI , G 

DIMENSION X(N) , FX(M) , TT(1584), XRI(1584) 

DIMENSION PH(4) , PP(4), PW(4), EQ(4) 

A1 = X(N-3) 

A2 = X( N-2) 

A3 = X(N-l) 

A4 = X(N) 
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DO 100 1=1, M 
INDEX = 0 
DO 200 L=1 ,4 
PH(L) = X( L+I NDEX) 

PP(L) = X(L+1+INDEX) 

PW(L) = X(L+2+INDEX) 

INDEX = INDEX + 2 

EQ(L) = (G *PH( L) * EXP(-AL0G(2.0)*((2.0*(TT(I)- PP(L)) / PW( 

1 L) )**2 . 0) ) )+ ( ( ( 1 .0-G)*PH(L) )/ (1.0+(((2.0*(TT(I)- PP(L))) / PW( 

2 L))**2.0))) 

200 CONTINUE 

B = A1 + (A2*TT( I) ) + (A3*(TT(I)**2)) + (A4*(TT( I)**3) ) 

50 CONTINUE 

FX( I) = EQ( 1) + EQ( 2) + EQ( 3) + EQ(4) - XRI(I) + B 
100 CONTINUE 
RETURN 
END 

SUBROUTINE AREA(X,M,N) 

COMMON TT, XRI, G 

DIMENSION TT( 1584) , XRI(1584), CP(1584) 

DIMENSION PH(4) , PP(4), PW(4), EQ(1584,4) 

DIMENSION IRAY(6) 

DIMENSION B(4) , X(16) 

DATA IRAY/6 * (-0)/ 

IRAY(4) = 0 

CALL SYSTEMC (115.IRAY) 

DO 110 1=1, M, 8 

PRINT6 , XRI( I) ,XRI(I+1) , XRI (1+2) ,XRI(I+3) , XRI (1+4) ,XRI( 1+5) , XRI( 1+ 

26) ,XRI(I+7) 

6 F ORMAT ( 1 H ,10X,8F10.2) 

110 CONTINUE 

DO 111 1=1, M, 8 

PRI NT6 ,TT( I) , TT( 1+1) ,TT(I+2) ,TT( 1+3) ,TT(I+4) ,TT(I+5) ,TT(I+6) ,TT(I+ 

27) 

111 CONTINUE 

B1 = X( N-3) 

B2 = X( N-2) 

B3 = X(N-l) 

B4 = X(N) 

DO 150 1=1, M 
INDEX = 0 
DO 200 L=1 , 4 
PH(L) = X(L+INDEX) 

PP( L) = X(L+1+INDEX) 

PW( L) = X(L+2+INDEX) 

INDEX = INDEX + 2 

EQ( I , L) = (G *PH( L) * EXP(-ALOG(2.0)*((2.0*(TT(I)- PP(L)) / PW( 

1 L) )**2. 0) ) )+ ( ( ( 1 . 0-G)*PH(L) )/ ( 1 .0+( ( ( 2 . 0*(TT( I)- PP(L))) / PW( 

2 L) )**2 . 0) ) ) 

200 CONTINUE 

CP(I) = B1 + B2 *TT( I) + B3 *(TT(I)**2) + B4 *(TT(I)**3) 

150 CONTINUE 


23 



COMMENT XRI= X-RAY INTENSITY; TAE=TOTAL AREA UNDER ENVELOPE 

TAE = 0.0 
A1 = 0.0 
A2 = 0.0 
A3 = 0.0 
A4 = 0.0 
A5 = 0.0 
PC = 0.0 
DO 225 1=1, M 
TAE = TAE + XRI(I) 

A1 = A1 + EQ( 1,1) 

A2 = A2 + EQ(I,2) 

A3 = A3 + EQ( I ,3) 

A4 = A4 + EQ( 1, 4) 

PC = PC + CP(I) 

225 CONTINUE 

COMMENT: GENERATE AREAS 

BASE = TT(M) - TT( I ) 

TAE = (TAE/M) * BASE 

A1 = (Al/M) * BASE 

A2 = (A2/M) * BASE 

A3 = (A3/M) * BASE 

A4 = (A4/M) * BASE 

PC = (PC/M) * BASE 

TOTAL = A1 + A2 + A3 + A4 +PC 

XTAL = ( A1 + A2 + A3 + A4) / TOTAL 


C 


PRINTIO, TAE, TOTAL 

10 FORMAT( 1H ,//, 1X.47HTOTAL AREA UNDER MEASURED INTENSITY ENVELOPE = 
1 , F14.4, 10X, I5H(T0TAL PEAKS = ,F14.4,1H)) 

PRINTI 1 ,A1 ,A2 , A3, A4, A5 

11 FORMAT( 1H0 , 28HAREAS OF DIFFRACTION PEAKS: ,5(F14.4,3X)) 

PRINT12 , PC 

12 FORMAT( 1H0, 34HAREA OF PARACRYSTALLINE SCATTER = ,F14.4) 

PRINTI 3, XTAL 

13 FORMAT( 1H0 , 23UCRYSTALLINE FRACTION = ,F14.4) 

RETURN 

END 


SUBROUTINE MARQ(M,N,X,C,EPS,NSIG,ITMAX,FX,F,WK,IERR) 

DIMENSION WK(5) ,X( 1) ,FX( 1) 

EXTERNAL F 

INDWKL = 5*N+2*M+(N+l)*N/2+l 
INDWK2 = INDWK1+1 
INDWK3 = INDWK2+M*N 
IOPT = 0 
DELTA = 0.0 

MAXFN = (ITMAX+1) * (N+l) 

CALL QXZ031(F, M,N, NSIG.EPS, DELTA, MAXFN, IOPT, C, X, WK(INDWKl) ,FX, 
*WK( INDWK2) ,M,WK(INDWK3),WK, INFER, IERR) 

ITMAX = WK(5) 

RETURN 

END 

SUBROUTINE QXZ03 1 (F ,M ,N,NSIG, EPS , DELTA, MAXFN, IOPT, PARM, X,SSQ, 


MARQ 2 
MARQ 90 
MARQ 91 
MARQ 92 
MARQ 96 
MARQ 97 
MARQ 98 
MARQ 99 
MARQ 100 
MARQ 101 
QXZNA270 
MARQ 103 
MARQ 104 
MARQ 105 
MARQ 106 
QXZNA273 


04 



n o 


c 

c 

c 


c 


c 


c 


c 


FX,XJAC,IXJAC,XJTJ,WK, INFER, IERR) 

THE SOLUTION X IS A STATIONARY POINT. 


DIMENSION 


REAL 

* 

* 

* 

* 

* 


EXTERNAL F 
DATA 
DATA 
DATA 

* 

ie 


X( N) , FX(M) , PARM( 5) ,XJAC( I) ,XJTJ( 1 ) ,WK(6) 

XJAC USED INTERNALLY IN PACKED FORM 
AL.CONS2 , DELTA , DNORM ,DSQ , EPS , 
ERL2,ERL2X,FX,FO,FOSQ,FOSQS4,G,HALF, 
XJTJ,HH,ONE,ONEP10,ONEP5,ONESFO,AX, 

PREC , REL , RHH , S IG, SQDIF , SSQ , S SQOLD, SUM .TEN, 
TENTH , X , XDIF , XHOLD , XJAC , UP , WK , ZERO , 

XDABS , RELCON, PARM , PO 1 , TWO , HUNTW , DELTA2 

SIG/14.4/ 

AX/0.1/ 

PO I , TENTH, HALF , ZERO , ONE, QNEP5 , TWO , 

TEN, HUNTW, ONEPIO/O. 01, 0.1, 0.5, 0.0, 

1. ,1.5,2. ,10. 0,1. 2E2.1.E10/ 

ERROR CHECKS 


IERR = 0 

IF (M.LE . 0 .OR.M.GT . IXJAC .OR. N.LE. 0 .OR . I0PT.LT.0 .OR . IOPT .GT. 2 
C .OR. PARM(l) .LE. 0.0) GO TO 305 

IMJC = IXJAC-M 
IF (I0PT.NE.2) GO TO 5 

IF ( PARM( 2) .LE .ONE .OR . PARM( 1 ) .LE . ZERO) GO TO 305 

MACHINE DEPENDENT CONSTANTS 


5 PREC = TEN** (-SIG-ONE) 

REL = TEN** (-SIG* HALF) 

RELCON = TEN** (-NSIG) 

WORK VECTOR IS CONCATENATION OF 
SCALED HESSIAN , GRADIENT , DELX, SCALE , 
XNEW,XBAD,F(X+DEL),F(X-DEL) 

IGRADI = ((N+l)*N)/2 
IGRADL = I GRAD 1+1 
IGRADU = IGRAD1+N 
IDELX1 = IGRADU 
IDELXL = IDELX1+1 
IDELXU = IDELX1+N 
ISCAL1 = IDELXU 
ISCALL = ISCAL1+1 
ISCALU = ISCAL1+N 
1XNEW1 = ISCALU 
IXNEWL =* IXNEW1+1 
IXBAD1 = IXNEW1+N 
IFPL1 =■ IXBAD1+N 
IF PL = IFPL1+1 
IFPU = IFPL1+M 
IFMLl =» IFPU 
IFML = IFML1+1 
IMJC = IXJAC - M 


INITIALIZE VARIABLES 


AL = PARM(l) 

C0NS2 = 0.0 

IF ( IOPT.EQ.O) GO TO 20 
IF (IOPT.EQ.l) GO TO 10 


ZXSSQL 3 
ZXSSQ150 
ZXSSQ170 
ZXSSQ171 
ZXSSQ172 
ZXSSQ173 
ZXSSQ174 
ZXSSQ175 
ZXSSQ176 
ZXSSQ177 
ZXSSQ178 
ZXSSQ179 
ZXSSQ180 
ZXSSQ181 
ZXSSQ182 
ZXSSQ183 
ZXSSQ184 
ZXSSQ185 
ZXSSQ186 
ZXSSQ187 
ZXSSQ188 
ZXSSQ189 
ZXSSQ190 
ZXSSQ191 
ZXSSQ192 
ZXSSQ193 
ZXSSQ194 
ZXSSQ195 
ZXSSQ196 
ZXSSQ197 
ZXSSQ198 
ZXSSQ199 
ZXSSQ200 
ZXSSQ201 
ZXSSQ202 
ZXSSQ203 
ZXSSQ204 
ZXSSQ205 
ZXSSQ206 
ZXSSQ207 
ZXSSQ208 
ZXSSQ209 
ZXSSQ210 
ZXSSQ21 1 
ZXSSQ212 
ZXSSQ213 
ZXSSQ214 
ZXSSQ215 
ZXSSQ216 
ZXSSQ217 
ZXSSQ218 
ZXSSQ219 
ZXSSQ220 
ZXSSQ221 
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AL = PARM(l) 

ZXSSQ222 


F0 =* PARM( 2) 

ZXSSQ223 


UP = PARM( 3 ) 

ZXSSQ224 


C0NS2 = PARM(4)*HALF 

ZXSSQ225 


GO TO 15 

ZXSSQ226 

10 

AL = P01 

ZXSSQ227 


FO = TWO 

ZXSSQ228 


UP - HUNTW 

ZXSSQ229 


C0NS2 = TENTH 

ZXSSQ230 

15 

ONESFO = ONE/FO 

ZXSSQ231 


FOSQ = FO*FO 

ZXSSQ232 


F0SQS4 = F0SQ**4 

ZXSSQ233 

20 

IEVAL = 0 

ZXSSQ234 


DELTA2 = DELTA*HALF 

ZXSSQ235 


ERL2 => 0NEP10 

ZXSSQ236 


IBAD = -99 

ZXSSQ237 


ISW = 1 

ZXSSQ238 


ITER = -1 

ZXSSQ239 


INFER = 0 

ZXSSQ240 


IERR = 0 

ZXSSQ241 


DO 25 J=IDELXL, IDELXU 

ZXSSQ242 


WK(J) = ZERO 

ZXSSQ243 

25 

CONTINUE 

ZXSSQ244 


GO TO 165 

ZXSSQ245 


MAIN LOOP 

ZXSSQ246 

30 

SSQOLD = SSQ 

ZXSSQ247 


CALCULATE JACOBIAN 

ZXSSQ248 


IF ( INFER .GT. O.OR. IJAC .GE .N.OR. IOPT.EQ. O.OR. ICOUNT.GT . 0) GO TO 

55 ZXSSQ249 


RANK ONE UPDATE TO JACOBIAN 

ZXSSQ250 


IJAC = IJAC+1 

ZXSSQ251 


DSQ = ZERO 

ZXSSQ252 


DO 35 J=IDELXL, IDELXU 

ZXSSQ253 


DSQ = DSQ + WK( J) * WK( J) 

ZXSSQ254 

35 

CONTINUE 

ZXSSQ255 


IF (DSQ. LE. ZERO) GO TO 55 

ZXSSQ256 


DO 50 1=1, M 

ZXSSQ257 


G = FX( I) - WK( IFML1 + I) 

ZXSSQ258 


K = I 

ZXSSQ259 


DO 40 J=IDELXL, IDELXU 

ZXSSQ260 


G = G + XJAC(K) * WK( J) 

ZXSSQ261 


K = K+IXJAC 

ZXSSQ262 

40 

CONTINUE 

ZXSSQ263 


G = G/DSQ 

ZXSSQ264 


K = I 

ZXSSQ265 


DO 45 J=IDELXL, IDELXU 

ZXSSQ266 


XJAC(K) = XJAC(K) - G * WK(J) 

ZXSSQ267 


K = K+IXJAC 

ZXSSQ268 

45 

CONTINUE 

ZXSSQ269 

50 

CONTINUE 

ZXSSQ270 


GO TO 80 

ZXSSQ27 1 


JACOBIAN BY INCREMENTING X 

ZXSSQ272 

55 

IJAC = 0 

ZXSSQ273 


K = -IMJC 

ZXSSQ274 


DO 75 J=1,N 

ZXSSQ275 


OR 



K = K+IMJC 
XDABS = ABS( X( J) ) 

HH = REL*(AMAX1( XDABS, AX) ) 
XHOLD = X(J) 

X(J) = X(J)+HH 
CALL F(X,M,N,WK( IFPL) ) 

IEVAL = IEVAL +1 
X(J) = XHOLD 
IF (ISW.EQ.l) GO TO 65 


C 

C 


CENTRAL DIFFERENCES 


60 


65 


70 

75 

80 


X(J) = XHOLD-HH 
CALL F(X,M,N, WK( IFML) ) 

IEVAL = IEVAL + 1 
X(J) =* XHOLD 
RHH = HALF/HH 
DO 60 I=IFPL, IFPU 
K = K+l 

XJAC(K) = ( WK( I) - WK( I 
CONTINUE 
GO TO 75 


+ M)) * RHH 


RHH = ONE/HH 
DO 70 1=1, M 
K = K+l 
XJAC(K) = 
CONTINUE 
CONTINUE 


FORWARD DIFFERENCES 


(WK( IFPL1 + I) - FX( I) ) * RHH 


85 


90 


ERL2X = ERL2 
ERL2 = ZERO 
K = -IMJC 

DO 90 J=IGRADL, IGRADU 
K = K+IMJC 
SUM = ZERO 
DO 85 1=1, M 
K = K+l 

SUM = SUM+XJAC(K)*FX( I) 
CONTINUE 
WK(J) = SUM 
ERL2 = ERL2+S UM*SUM 
CONTINUE 

ERL2 = SQRT(ERL2) 

IF ( IJAC.GT.O) GO TO 95 
IF ( ERL2 . LE . DELTA2 ) INFER = 

IF (ERL2.LE.CONS2) ISW = 2 


95 


L = 0 

IS = -IXJAC 
DO 110 1=1, N 

IS = I S+ IXJAC 
JS = -IXJAC 
DO 105 J=1,I 


ZXSSQ276 
ZXSSQ277 
ZXSSQ278 
ZXSSQ279 
ZXSSQ280 
ZXSSQ281 
ZXSSQ282 
ZXSSQ283 
ZXSSQ284 
ZXSSQ285 
ZXSSQ286 
ZXSSQ287 
ZXSSQ288 
ZXSSQ289 
ZXSSQ290 
ZXSSQ291 
ZXSSQ292 
ZXSSQ293 
ZXSSQ294 
ZXSSQ295 
ZXSSQ296 
ZXSSQ297 
ZXSSQ298 
ZXSSQ299 
ZXSSQ300 
ZXSSQ301 
ZXSSQ302 
ZXSSQ303 
ZXSSQ304 
ZXSSQ305 
ZXSSQ306 
ZXSSQ307 
ZXSSQ308 
ZXSSQ309 
ZXSSQ310 
ZXSSQ311 
ZXSSQ312 
ZXSSQ313 
ZXSSQ314 
ZXSSQ315 
ZXSSQ316 
ZXSSQ317 

CONVERGENCE TEST FOR NORM OF GRADIENTZXS SQ3 18 

ZXSSQ319 

INFER+4 ZXSSQ320 

ZXSSQ321 
ZXSSQ322 
ZXSSQ323 
ZXSSQ324 
ZXSSQ325 
ZXSSQ326 
ZXSSQ327 
ZXSSQ328 
ZXSSQ329 


CALCULATE GRADIENT 


CALCULATE THE LOWER SUPER TRIANGE OF 
JACOBIAN (TRANSPOSED) * JACOBIAN 


27 




JS = JS+IXJAC 

ZXSSQ330 


L = L+l 

ZXSSQ331 


SUM = ZERO 

ZXSSQ332 


DO 100 K=1,M 

ZXSSQ333 


LI = IS+K 

ZXSSQ334 


LJ = JS+K 

ZXSSQ335 


SUM - SUM+XJAC( LI) *X JAC( LJ) 

ZXSSQ336 

100 

CONTINUE 

ZXSSQ337 


XJTJ(L) = SUM 

ZXSSQ338 

105 

CONTINUE 

ZXSSQ339 

no 

CONTINUE 

ZXSSQ340 


CONVERGENCE CHECKS 

ZXSSQ341 


IF ( INFER. GT.O) GO TO 315 

ZXSSQ342 


IF( IEVAL .GE. MAXFN) GO TO 295 

ZXSSQ343 


COMPUTE SCALING VECTOR 

ZXSSQ344 


IF (IOPT.EQ.O) GO TO 120 

ZXSSQ345 


K = 0 

ZXSSQ346 


DO 115 J=1,N 

ZXSSQ347 


K = K+J 

ZXSSQ348 


WK( ISCAL1 + J) = XJTJ(K) 

ZXSSQ349 

115 

CONTINUE 

ZXSSQ350 


GO TO 135 

ZXSSQ351 


COMPUTE SCALING VECTOR AND NORM 

ZXSSQ352 

120 

DNORM = ZERO 

ZXSSQ353 


K = 0 

ZXSSQ354 


DO 125 J=1,N 

ZXSSQ355 


K = K+J 

ZXSSQ356 


WK( ISCAL1 + J) = SQRT(XJTJ(K)) 

ZXSSQ357 


DNORM = DNORM+XJTJ( K)*XJTJ(K) 

ZXSSQ358 

125 

CONTINUE 

ZXSSQ359 


DNORM = ONE/ SQRT( DNORM) 

ZXSSQ360 


NORMALIZE SCALING VECTOR 

ZXSSQ361 


DO 130 J=ISCALL, ISCALU 

ZXSSQ362 


WK( J) = WK( J) * DNORM * ERL2 

ZXSSQ363 

130 

CONTINUE 

ZXSSQ364 


ADD L-M FACTOR TO DIAGONAL 

ZXSSQ365 

135 

ICOUNT = 0 

ZXSSQ366 

140 

K = 0 

ZXSSQ367 


DO 150 1=1, N 

ZXSSQ368 


DO 145 J=1,I 

ZXSSQ369 


K = K+l 

ZXSSQ370 


WK( K) = XJTJ(K) 

ZXSSQ37 1 

145 

CONTINUE 

ZXSSQ372 

147 

WK(K) = WK( K) + WK(ISCAL1 + I) * AL 

ZXSSQ373 

148 

WK( IDELX1 + I) = WK( IGRAD1 + I) 

ZXSSQ374 

150 

CONTINUE 

ZXSSQ375 


CHOLESKY DECOMPOSITION 

ZXSSQ376 

155 

CALL QXZ032 (WK, 1 , N, WK( IDELXL) ,N,0,G,XH0LD,IERR) 

QXZNA276 


IF( IERR.EQ. 0) GO TO 160 

ZXSSQ378 


IERR = 0 

ZXSSQ379 


IF ( IJAC.GT. 0) GO TO 55 

ZXSSQ380 


IF (IBAD.LE.O) GO TO 240 

ZXSSQ381 


IF (IBAD.GE.2) GO TO 310 

ZXSSQ382 


GO TO 190 

ZXSSQ383 
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160 IF (IBAD.NE.-99) IBAD = 0 

C CALCULATE SUM OF SQUARES 

165 DO 170 J=1,N 

WK( IXNEW1 + J) = X(J) - WK( IDELX1 + J) 

170 CONTINUE 

CALL F(WK( IXNEWL) ,M,N, WK( I FPL) ) 

IEVAL = IEVAL + 1 

SSQ = ZERO 

DO 175 I=IFPL, IFPU 

SSQ = SSQ + WK( I) * WK( I) 

175 CONTINUE 

177 IF (ITER.GE.O) CO TO 185 

C SSQ FOR INITIAL ESTIMATES OF X 

ITER = 0 
SSQOLD = SSQ 
DO 180 1=1, M 

FX( I) = WK( IFPL1 + I) 

180 CONTINUE 
GO TO 55 

185 IF ( IOPT.EQ. 0) GO TO 215 

C CHECK DESCENT PROPERTY 

IF (SSQ. LE. SSQOLD) GO TO 205 

C INCREASE PARAMETER AND TRY AGAIN 

190 ICOUNT = ICOUNT+1 
AL = AL*FOSQ 
IF ( IJAC.EQ.O) GO TO 195 
IF (IC0UNT.GE.4.0R.AL.GT.UP) GO TO 200 
195 IF (AL.LE.UP) GO TO 140 
IF ( IBAD.EQ. 1) GO TO 310 
GO TO 300 

200 AL = AL/F0SQS4 
GO TO 55 

C ADJUST MARQUARDT PARAMETER 

205 IF ( ICOUNT. EQ.O) AL = AL/FO 
IF (ERL2X.LE. ZERO) GO TO 210 
G = ERL2/ERL2X 

IF ( ERL2.LT. ERL2X) AL = AL*AMAX1(ONESFO,G) 

IF (ERL2.GT.ERI.2X) AL = AL*AMIN1(F0,G) 

210 AL = A MAX 1 ( AL , PREC) 

C ONE ITERATION CYCLE COMPLETED 

215 ITER = ITER+1 
DO 220 J=1,N 

X(J) = WK( IXNEW1 + J) 

220 CONTINUE 

DO 225 1=1, M 

WK( IFML1 + I) = FX( I) 

FX( I) = WK( IFPL1 + I) 

225 CONTINUE 

C RELATIVE CONVERGENCE TEST FOR X 

IF (ICOUNT. GT.O. OR. IJAC.GT.O) GO TO 30 
DO 230 J=1,N 

XDIF = ABS( WK( IDELX1 + J)) / AMAX1(ABS(X( J) ) ,AX) 

IF (XDIF.GT.RELCON) GO TO 235 
230 CONTINUE 


ZXSSQ384 

ZXSSQ385 

ZXSSQ386 

ZXSSQ387 

ZXSSQ388 

ZXSSQ389 

ZXSSQ390 

ZXSSQ391 

ZXSSQ392 

ZXSSQ393 

ZXSSQ394 

ZXSSQ395 

ZXSSQ396 

ZXSSQ397 

ZXSSQ398 

ZXSSQ399 

ZXSSQ400 

ZXSSQ401 

ZXSSQ402 

ZXSSQ403 

ZXSSQ404 

ZXSSQ405 

ZXSSQ406 

ZXSSQ407 

ZXSSQ408 

ZXSSQ409 

ZXSSQ410 

ZXSSQ411 

ZXSSQ412 

ZXSSQ413 

ZXSSQ414 

ZXSSQ415 

ZXSSQ416 

ZXSSQ417 

ZXSSQ418 

ZXSSQ419 

ZXSSQ420 

ZXSSQ421 

ZXSSQ422 

ZXSSQ423 

ZXSSQ424 

ZXSSQ425 

ZXSSQ426 

ZXSSQ427 

ZXSSQ428 

ZXSSQ429 

ZXSSQ430 

ZXSSQ431 

ZXSSQ432 

ZXSSQ433 

ZXSSQ434 

ZXSSQ435 

ZXSSQ436 

ZXSSQ437 





n n 


INFER = INFER + 2 

C RELATIVE CONVERGENCE TEST FOR SSQ 

235 SQDIF = ABS(SSQ-SSQOLD)/AMAXl( SSQOLD, AX) 

IF (SQDIF .LE. EPS) INFER - INFER+1 
IF( IBAD .EQ. -99) GO TO 30 

IF( INFER .NE. 0 .AND. SSQ .GT. (10.0 * EPS)) GO TO 310 
GO TO 30 

C SINGULAR DECOMPOSITION 

240 IF (IBAD) 255,245,265 

CHECK TO SEE IF CURRENT 
ITERATE HAS CYCLED BACK TO 

C THE LAST SINGULAR POINT 

245 DO 250 J=1,N 

XHOLD = WK( IXBAD1 + J) 

IF (ABS(X(J)-XH0LD).GT.RELC0N*AMAX1(AX,ABS(XH0LD))) GO TO 255 
250 CONTINUE 
GO TO 310 

C UPDATE THE BAD X VALUES 

255 DO 260 J=1,N 

WK( IXBADI + J) = X(J) 

260 CONTINUE 
IBAD = 1 

C INCREASE DIAGONAL OF HESSIAN 

265 IF ( IOPT.NE. 0) GO TO 280 
K = 0 

DO 275 1=1, N 
DO 270 J=1,I 
K = K+l 

WK(K) = XJTJ(K) 

270 CONTINUE 

WK(K) = 0NEP5 * (XJTJ(K) + AL * ERL2 * WK(ISCAL1 + I)) + REL 
275 CONTINUE 
IBAD = 2 
GO TO 155 

C REPLACE ZEROES ON HESSIAN DIAGONAL 


280 


285 


290 

295 

300 

305 

310 


C 

315 


IZERO = 0 

DO 285 J=ISCALL, ISCALU 

IF( WK( J) . GT . ZERO) GO TO 285 
IZERO = IZERO+1 
WK( J) = ONE 
CONTINUE 

IF (IZERO. LT.N) GO TO 140 

TERMINAL ERROR 

IERR = IERR + 1 
IERR = IERR + 1 
IERR = IERR + 1 
IERR = IERR + 1 
CONTINUE 


IERR = IERR + 1 

IF (IERR . EQ . 2) GO TO 9005 

OUTPUT ERL2 , IEVAL,NSIG, AL, AND ITER 

G = SIG 
DO 320 J=1 , N 

XHOLD = AB S ( WK( IDELX1 + J)) 


ZXSSQ438 

ZXSSQ439 

ZXSSQ440 

ZXSSQ441 

ZXSSQ442 

ZXSSQ443 

ZXSSQ444 

ZXSSQ445 

ZXSSQ446 

ZXSSQ447 

ZXSSQ448 

ZXSSQ449 

ZXSSQ450 

ZXSSQ451 

ZXSSQ452 

ZXSSQ453 

ZXSSQ454 

ZXSSQ455 

ZXSSQ456 

ZXSSQ457 

ZXSSQ458 

ZXSSQ459 

ZXSSQ460 

ZXSSQ461 

ZXSSQ462 

ZXSSQ463 

ZXSSQ464 

ZXSSQ465 

ZXSSQ466 

ZXSSQ467 

ZXSSQ468 

ZXSSQ469 

ZXSSQ470 

ZXSSQ471 

ZXSSQ472 

ZXSSQ473 

ZXSSQ474 

ZXSSQ475 

ZXSSQ476 

ZXSSQ477 

ZXSSQ478 

ZXSSQ479 

ZXSSQ480 

ZXSSQ481 

ZXSSQ482 

ZXSSQ483 

ZXSSQ484 

ZXSSQ485 

ZXSSQ486 

ZXSSQ487 

ZXSSQ488 

ZXSSQ489 

ZXSSQ490 

ZXSSQ491 


30 



n o 



IF (XHOLD.LE.ZERO) GO TO 320 

ZXSSQ492 


G = AMIN1(G, -ALOG10(XH0LD)+ALOG10( AMAX1 ( AX,ABS(X( J) ) ) ) ) 

ZXSSQ493 

320 

CONTINUE 

ZXSSQ494 


IF(N.GT. 2) GO TO 330 

ZXSSQ495 


DO 325 J = 1, N 

ZXSSQ496 

325 

WK( J+5) => WK( J+IGRAD1) 

ZXSSQ497 

330 

WK( I) = ERL2+ERL2 

ZXSSQ498 


WK( 2) = IEVAL 

ZXSSQ499 


SSQ = SSQOLD 

ZXSSQ500 


WK( 3) = G 

ZXSSQ501 


WK(4) = AL 

ZXSSQ502 


WK(5) = ITER 

ZXSSQ503 

9005 

RETURN 

ZXSSQ504 


END 

ZXSSQ505 


SUBROUTINE QXZ032 (A,M,N, B , IB , IDGT.Dl , D2 , IER) 

QXZNA2 7 9 
LEQT1P 3 

; LANGUAGE - FORTRAN 

LEQT1P39 

LEQT1P42 


DIMENSION A( 1) ,B(IB, I) 

LEQT1P43 


INITIALIZE IER 

LEQT1P44 


IER = 0 

LEQT1P45 


DECOMPOSE A 

LEQT1P46 


CALL QXZ033 (A,A,N,DL,D2,IER) 

QXZNA282 


IF (IER.NE.O) GO TO 9005 

LEQT1P48 


PERFORM ELIMINATION 

LEQT1P49 


DO 5 I = 1,M 

LEQT1P50 


CALL QXZ034 (A,B(1,I) ,N,B(1,I)) 

QXZNA283 

5 

CONTINUE 

LEQT1P52 

9005 

RETURN 

LEQT1P53 


END 

LEQT1P54 


SUBROUTINE QXZ033 (A,UL,N,D1,D2,IER) 

QXZNA286 
LUDECP 3 

: PRECISION - SINGLE 

LUDECP25 


DIMENSION A(1),UL(1) 

LUDECP31 


DATA ZERO, ONE, FOUR, SIXTN, SIXTH/O. 0,1. ,4. ,16. ,.0625/ 

LUDECP32 


D1=0NE 

LUDECP33 


D2=ZER0 

LUDECP34 


RN = ONE/ (N*SIXTN) 

LUDECP35 


IP = I 

LUDECP 3 6 


IER=0 

LUDECP37 


DO 45 I “ 1,N 

LUDECP 3 8 


IQ = IP 

LUDECP39 


IR = 1 

LUDECP40 


DO 40 J = 1,1 

LUDECP4 1 


X = A( IP) 

LUDECP42 


IF (J .EQ. 1) GO TO 10 

LUDECP43 


DO 5 K=*IQ,IP1 

LUDECP44 


X = X-UL( K)*UL( IR) 

LUDECP45 


IR = IR+1 

LUDECP46 

5 

CONTINUE 

LUDECP47 

10 

IF ( I.NE.J) GO TO 30 

LUDECP48 


D1 = DI*X 

LUDECP49 


IF (A(IP)+X*RN .LE. A( IP) ) GO TO 50 

LUDECP50 

15 

IF (ABS(Dl) .LE. ONE) GO TO 20 

LUDECP51 




D1 = D1 * SIXTH 

LUDECP52 


D2 * D2 + FOUR 

LUDECP53 


GO TO 15 

LUDECP54 

20 

IF (ABS(Dl) .GE. SIXTH) GO TO 25 

LUDECP55 


D1 = Dl * SIXTN 

LUDECP56 


D2 = D2 - FOUR 

LUDECP57 


GO TO 20 

LUDECP58 

25 

UL(IP) - ONE/ SQRT(X) 

LUDECP59 


GO TO 35 

LUDECP60 

30 

UL(IP) = X * UL(IR) 

LUDECP61 

35 

IP1 = IP 

LUDECP62 


IP = IP+1 

LUDECP63 


IR = IR+1 

LUDECP64 

40 

CONTINUE 

LUDECP65 

45 

CONTINUE 

LUDECP66 


GO TO 9005 

LUDECP67 

50 

IER = 129 

LUDECP68 

9000 

CONTINUE 

LUDECP69 

9005 

RETURN 

LUDECP7 0 


END 

LUDECP7 1 


SUBROUTINE QXZ034 (A,B,N,X) 

QXZNA290 



LUELMP 3 



LUELMP25 


DIMENSION A( 1) ,B( 1) ,X( 1) 

LUELMP26 


DATA ZERO/O. / 

LUELMP 27 


SOLUTION OF LY = B 

LUELMP28 


IP=1 

LUELMP 2 9 


IW = 0 

LUELMP 30 


DO 15 1=1, N 

LUELMP3 1 


T=B( I) 

LUELMP32 


IM1 = 1-1 

LUELMP 3 3 


IF (IW .EQ. 0) GO TO 9 

LUELMP34 


IP=IP+IW-1 

LUELMP 3 5 


DO 5 K= IW , IM1 

LUELMP36 


T = T-A( IP)*X(K) 

LUELMP 3 7 
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GO TO 10 
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LUELMP54 


KK = N 

LUELMP 5 5 
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X( II)=T*A( IS) 
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CONTINUE 
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